clear variables
datapath    = 'Output\';
filename    = 'Chains.xlsx';
name        = 'replication';
cost        = 1;
aux.n_m     = 10^3;
aux.n       = 10^3;

% loop over specifications
for shock_negative = 0:5
        
    path='parameters';
    
    if shock_negative == 0 
        run_spec = '_OJS';
    elseif shock_negative == 1 
        run_spec = '_27';
    elseif shock_negative == 2 
        run_spec = '_10';
    elseif shock_negative == 3 
        run_spec = '_38';
    elseif shock_negative == 4
        run_spec = '_17';
    elseif shock_negative == 5
        run_spec = '_22';
    end
    if shock_negative == 5
        sheet_name = 'Table 1, Table 3';
    else
        sheet_name = ['Parameters',run_spec];
    end

    % load parameters
    parampath = [path,run_spec];
    load([parampath, '\param_ini_', name, '.mat']);
    load([parampath, '\param_end_', name, '.mat']);
    L=20./(1-u_fun(param_ini));% calculate implied labor force
        
    %% Save parameters
    varNames = { 'alpha', 'r', 'sigma', ...
        'lambda', 's', 'exo_sep', 'beta', 'omega_bar', ...
        'omega_0', 'L', 'c', 'K', 'mean_x'};

    param=param_ini;

    T = table(param_ini.alpha, param_ini.r, param_ini.sigma, ...
    param_ini.lambda, param_ini.s, param_ini.zeta, param_ini.beta, log(param_end.omega_0/param_ini.omega_0)./log(1-param_ini.shock), ...
    param_ini.omega_0, L, param_ini.c, param_ini.K, param_ini.mean_x*L, ...
                 'VariableNames',varNames);
    writetable(T,[datapath,filename],'Sheet',sheet_name,'Range','A1')

    % Calculate moments in Table 1, 3 (and Table 4 for alternative
    % specifications)
    m = linspace(param.m_l,param.m_u,10^3)';
    q = q_fun(m,param);
    G = (q-q(1))./(q(end)-q(1));
    G_r = (q_fun(param.m_h,param)-q(1))./(q(end)-q(1));

    dlnu_dlnAPL = log(u_fun(param_ini)./u_fun(param_end))./log(mean_fun(aux.n,'m',param_ini)./mean_fun(aux.n,'m',param_end));
    
    dlnlam_dlnAPL = log(param_ini.lambda./param_end.lambda)./log(mean_fun(aux.n,'m',param_ini)./mean_fun(aux.n,'m',param_end));
    
    dlnEU_dlnAPL = log(u_fun(param_ini).*param_ini.lambda/(1-u_fun(param_ini)) /...
        (u_fun(param_end).*param_end.lambda/(1-u_fun(param_end))))./log(mean_fun(aux.n,'m',param_ini)./mean_fun(aux.n,'m',param_end));
    
    m_ini       = linspace(param_ini.m_l,param_ini.m_u,10^3)';
    m_end       = linspace(param_end.m_l,param_end.m_u,10^3)';

    q_ini       = q_fun(m_ini,param_ini);
    q_p_ini     = q_ini(2:end)-q_ini(1:end-1);
    q_p_ini     = q_p_ini./sum(q_p_ini);
    E_lnw_ini   = q_p_ini'*log(wage_fun( m_ini(2:end), param_ini));
    
    q_end       = q_fun( m_end, param_end);
    q_p_end     = q_end(2:end) - q_end(1:end-1);
    q_p_end     = q_p_end./sum(q_p_end);
    E_lnw_end   = q_p_end'*log(wage_fun( m_end(2:end), param_end));
    
    dElnw_du = (E_lnw_ini - E_lnw_end)./(u_fun(param_ini) - u_fun(param_end));
    
    
    dlnEE_dlnAPL = log( jtj_num( delta_fun( m_ini, param_ini), q_fun( m_ini, param_ini))/ ...
        jtj_num( delta_fun( m_end, param_end), q_fun( m_end, param_end)))./log( mean_fun(aux.n,'m',param_ini)./mean_fun(aux.n,'m',param_end));
    
    chain_ini =  (jtj_num( delta_fun( m_ini, param_ini), q_fun( m_ini, param_ini)) *(1 - u_fun(param_ini)) + param_ini.lambda*u_fun(param_ini)) ...
        ./(param_ini.s*param_ini.lambda*(q_fun(param_ini.m_h,param_ini)-q_fun(param_ini.m_l,param_ini))./(q_fun(param_ini.m_u,param_ini)-q_fun(param_ini.m_l,param_ini)).*(1 - u_fun(param_ini)) + param_ini.lambda*u_fun(param_ini));
    
    chain_end =  (jtj_num( delta_fun( m_end, param_end), q_fun( m_end, param_end)) *(1 - u_fun(param_end)) + param_end.lambda*u_fun(param_end)) ...
        ./(param_end.s*param_end.lambda*(q_fun(param_end.m_h,param_end)-q_fun(param_end.m_l,param_end))./(q_fun(param_end.m_u,param_end)-q_fun(param_end.m_l,param_end)).*(1 - u_fun(param_end)) + param_end.lambda*u_fun(param_end));
    
    dlnChain_dlnAPL = log(chain_ini./chain_end)./log(mean_fun(aux.n,'m',param_ini)./mean_fun(aux.n,'m',param_end));
    
    hires_rate_ini =  (jtj_num( delta_fun( m_ini, param_ini), q_fun( m_ini, param_ini)) *(1 - u_fun(param_ini)) + param_ini.lambda*u_fun(param_ini))./(1 - u_fun(param_ini));
    
    hires_rate_end =  (jtj_num( delta_fun( m_end, param_end), q_fun( m_end, param_end)) *(1 - u_fun(param_end)) + param_end.lambda*u_fun(param_end))./(1 - u_fun(param_end));
    
    dlnHires_rate_dlnAPL = log(hires_rate_ini./hires_rate_end)./log(mean_fun(aux.n,'m',param_ini)./mean_fun(aux.n,'m',param_end));
    
    T = table( chain_ini, dlnlam_dlnAPL, dlnEU_dlnAPL, dlnu_dlnAPL, dlnEE_dlnAPL,dlnHires_rate_dlnAPL,dlnChain_dlnAPL,dElnw_du);
    writetable(T,[datapath,filename],'Sheet',sheet_name,'Range','A4')

    % calculate moments in calibration
    if  param.K>0
    
        T_y = KFE_firms(param,aux);               
        ln_m=linspace(log(param.m_l),log(param.m_u),aux.n_m)';
        T_yy=T_y.*repmat((exp(ln_m)>param.m_h).*(exp(ln_m)<param.m_e),1,aux.n_m);
        T_yy=T_yy.*repmat((exp(ln_m')>param.m_h).*(exp(ln_m')<param.m_e),aux.n_m,1);
    
        q = q_fun(exp(ln_m),param);             
        q_p= q(2:end)-q(1:end-1);
        h1=[q_p;0]./sum(q_p);
        Tr=expm(T_yy);
        Tr=Tr.*repmat((exp(ln_m)>param.m_h).*(exp(ln_m)<param.m_e),1,aux.n_m);
        inac=sum(Tr^12*h1);
    end
      
    u=u_fun(param); 
    
    n=10^3;
    m           = linspace(param.m_l,param.m_u,n)'; 
    q           = q_fun(m,param);%-q_fun(param.m_l,param))./(q_fun(param.m_u,param)-q_fun(param.m_l,param));
    q_p         = q(2:end)-q(1:end-1);
    q_p         = q_p./sum(q_p);
    job_to_job  = jtj_num( delta_fun(m,param), q);
    E_w         = q_p'*wage_fun( m(2:end), param);
    
    delta = delta_fun(m,param);
    delta_p = -(delta(2:end)-delta(1:end-1));
    
    move = q_p'*cumsum(delta_p,'reverse');
    ma = (q_p./wage_fun(m(1:end-1),param))'*cumsum(delta_p.*wage_fun(m(1:end-1),param),'reverse');
    gain = ma./move-1;
    mm = repmat(m(1:end-1),1,length(m(1:end-1)));
    log_gain = sum(sum((delta_p.*q_p'.* tril(log(wage_fun(mm,param)./wage_fun(mm',param))))))./sum(sum(tril((delta_p.*q_p'))));
    c_Ew=param.c./E_w;
    
    ela_rent = q_p'*(log((wage_fun(m(2:end),param)./wage_fun(m(1:end-1),param)))./log((m(2:end)./m(1:end-1))));
    var_lnw = q_p'*(log(wage_fun(m(1:end-1),param))-q_p'*log(wage_fun(m(1:end-1),param))).^2;
    avrg_firm_size = L*(1-u_fun(param_ini));

    if param.K>0
        T = table( job_to_job, u, c_Ew, log_gain, inac, ela_rent, avrg_firm_size);
    else
        T = table( job_to_job, u, c_Ew, log_gain, ela_rent, avrg_firm_size);
    end     
    
    writetable(T,[datapath,filename],'Sheet',sheet_name,'Range','A7')    
    
    clear param*
end



for z = 1:2
    if z == 1
        run_spec = '_22';
    else
        run_spec = '_OJS';
    end
    clear param*
    parampath = [path,run_spec];

    load([parampath, '\param_ini_', name, '.mat']);
    load([parampath, '\param_end_', name, '.mat']);
    param=param_ini;

    if param.K>0
        

        T_y = KFE_firms(param,aux); 
        ln_m=linspace(log(param.m_l),log(param.m_u),aux.n_m)';
        T_yy=T_y.*repmat((exp(ln_m)>param.m_h).*(exp(ln_m)<param.m_e),1,aux.n_m);
        T_yy=T_yy.*repmat((exp(ln_m')>param.m_h).*(exp(ln_m')<param.m_e),aux.n_m,1);

        q   = q_fun(exp(ln_m),param);             
        q_p = q(2:end)-q(1:end-1);
        h1  = [q_p;0]./sum(q_p);
        Tr  = expm(T_yy);
        Tr  = Tr.*repmat((exp(ln_m)>param.m_h).*(exp(ln_m)<param.m_e),1,aux.n_m);

        T       = expm(10^3.*T_y);
        g_ini   = T*(1/aux.n_m.*ones(aux.n_m,1));
        g_ini   = T*g_ini;

        dt    = 1/30;
        tic
        Tr_dt = expm(T_yy.*dt);
        toc
        Tr_dt = Tr_dt.*repmat((exp(linspace(log(param.m_l),log(param.m_u),aux.n_m)')>param.m_h).*(exp(linspace(log(param.m_l),log(param.m_u),aux.n_m)')<param.m_e),1,aux.n_m);
        N = 48;
        inac = nan(N,1);
        inac_2 = nan(N,1);
        del_tmp = nan(N,1);
        delta_sum=nan(length(ln_m),N);
        h=nan(length(ln_m),N);

        for i=1:48
            g_ini=Tr*g_ini;
            if i==1
                h(:,i)=Tr*h1;
            else
                h(:,i)=Tr*h(:,i-1);
            end
            inac(i)      = sum(h(:,i));
            inac_2(i)      = sum(g_ini);
            if i==1
                delta_sum(:,i)=0;
            else
                delta_sum(:,i)=Tr*delta_sum(:,i-1);
            end

            for j=1:1:1/dt
                if i==1
                    delta_sum(:,i)=delta_sum(:,i)+(Tr_dt^(1/dt-j))*((max(param.zeta+delta_fun(exp(ln_m),param),0)).*(Tr_dt^(j)*h1)).*dt;
                else
                    delta_sum(:,i)=delta_sum(:,i)+(Tr_dt^(1/dt-j))*((max(param.zeta+delta_fun(exp(ln_m),param),0)).*(Tr_dt^(j)*h(:,i-1))).*dt;
                end
            end

            del_tmp(i)=sum(delta_sum(:,i))./inac(i);

        end
        Months = (1:48)';

        T = expm(T_yy);    
        for i = 1:aux.n_m
            g_tmp = zeros(size(g_ini));
            g_tmp(i) = 1;
            g_new = T*g_tmp;    

        end

        if param.K>0
            Cumulative_separations_C=del_tmp;
            inaction_emp_C=inac;
            inaction_firms_C=inac_2;
            T = table(Months,inaction_emp_C,Cumulative_separations_C);
            writetable(T,[datapath,filename],'Sheet','Figure 6','Range','A1')
        end
    end        

    param=param_ini;
    m = linspace(param.m_l,param.m_u,10^3)';
    J = J_fun(m,param);
    J_p = J_p_fun(m,param);
    J_pp = J_pp_fun(m,param);
    q   = q_fun(m,param);
    q   = q./q(end);
    q_p = q_p_fun(m,param);
    delta_p        = delta_fun(m(2:end),param)-delta_fun(m(1:end-1),param);
    g_p        = q(2:end)-q(1:end-1);
    g_p        = g_p./sum(g_p);
    job_to_job = g_p'*delta_fun(m(2:end),param);
    hires      = job_to_job.*(1-u_fun(param)) + param.lambda.*u_fun(param);
    hire_per_firm=sum(-delta_p.*q(1:end-1)./(q(end)-q(1)));

    if param.K>0
        rep_cost=rep_cost_fun( m, param);
    end

    delta = delta_fun(m,param);
    delta_p = delta_p_fun(m,param);
    eta = eta_fun(m,param);
    wage_m = wage_fun(m,param);
    T = table(m, J, J_p, J_pp, q, q_p, delta, delta_p, eta, wage_m);
    if param.K>0
        T = table(m, J, J_p, J_pp, q, q_p, delta, delta_p, eta, wage_m,rep_cost);
        writetable(T,[datapath,filename],'Sheet','Figure 4','Range','A1')
    end

    %% JC/BC

    n_j=30;
    lambda=linspace(0.1,0.5,n_j)';
    u_ED=nan(n_j,1);
    u_BC=nan(n_j,1);
    u_BC2=nan(n_j,1);
    u_ED2=nan(n_j,1);
    theta0=0.1;

    for j=1:n_j
        param_ini.lambda=lambda(j);
        param_end.lambda=lambda(j);

        param_ini=model(param_ini,aux);                
        param_end=model(param_end,aux);                

        u_BC(j)=u_fun(param_ini);
        u_BC2(j)=u_fun(param_end);

        mx  = mean_fun( aux.n, 'x', param_ini);
        mx2 = mean_fun( aux.n, 'x', param_end);

        e=param.mean_x/mx.*(1-u_fun(param_ini));
        u_ED(j)=1-e; 

        e=param.mean_x/mx2.*(1-param.shock).^(1/(1-param.alpha)).*(1-u_fun(param_end));
        u_ED2(j)=1-e; 

    end

    T = table(lambda,u_BC,u_BC2, u_ED, u_ED2);
    if param.K>0
            writetable(T,[datapath,filename],'Sheet','Figure 7','Range','A1')
    else
            writetable(T,[datapath,filename],'Sheet','Figure 7 (C=0 data)','Range','A1')
    end
end
